home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Fritz: All Fritz
/
All Fritz.zip
/
All Fritz
/
FILES
/
PROGMISC
/
FORTRAN1.LZH
/
INTRPL.FOR
< prev
next >
Wrap
Text File
|
1988-02-08
|
4KB
|
166 lines
SUBROUTINE INTRPL ( L, X, Y, N, U, V, IERR )
C
C INTRPL INTERPOLATES A FUNCTION Y(X) FROM A SET OF VALUES MONOTONICALY
C INCREASING IN X.
C
C INPUTS :
C L = NUMBER OF DATA POINTS (>= 2)
C X = ARRAY OF X VALUES IN INCREASING ORDER
C Y = ARRAY OF Y VALUES
C N = NUMBER OF POINTS TO BE INTERPOLATED
C U = ARRAY OF X LOCATIONS FOR OUTPUT (N VALUES)
C
C OUTPUTS :
C V = ARRAY OF INTERPOLATED Y VALUES
C IERR = ERROR INDICAT0R...
C 0 = NO ERROR
C 1 = L TOO SMALL
C 2 = N LESS OR EQUAL TO 0
C 3 = X VALUES NOT MONOTONICALY INCREASING
C
DIMENSION X(L), Y(L), U(N), V(N)
REAL M1,M2,M3,M4,M5
EQUIVALENCE (X2,A1,M1),(X5,A5,M5),(Q1,T3)
EQUIVALENCE (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3)
C
IERR = 0
LO = L
LM1 = LO - 1
LM2 = LM1 - 1
LP1 = LO + 1
NO = N
IF ( LM2 .LT. 0 ) THEN
IERR = 1
RETURN
ENDIF
IF ( NO .LE. 0 ) THEN
IERR = 2
RETURN
ENDIF
DO 11 I = 2, LO
IF(X(I-1) .GE. X(I)) THEN
IERR = 3
RETURN
ENDIF
11 CONTINUE
IPV = 0
DO 80 K = 1, NO
DX = U(K)
IF (LM2 .EQ. 0) THEN
I = 2
ELSE IF (DX .GE. X(LO)) THEN
I = LP1
ELSE IF (DX .LT. X(1)) THEN
I = 1
ELSE
IMN = 2
IMX = LO
21 I = (IMN + IMX)/2
IF (DX .LT. X(I)) THEN
IMX = I
ELSE
IMN = I + 1
ENDIF
IF (IMX .GT. IMN) GO TO 21
I = IMX
ENDIF
IF (I .EQ. IPV) GO TO 70
IPV = I
J = I
IF (J .EQ. 1)J=2
IF (J .EQ. LP1)J=LO
X3 = X(J-1)
Y3 = Y(J-1)
X4 = X(J)
Y4 = Y(J)
A3 = X4 - X3
M3 = (Y4 - Y3)/A3
IF (LM2 .EQ. 0) THEN
M2 = M3
M4 = M3
GO TO 45
ENDIF
IF (J .NE. 2) THEN
X2 = X(J-2)
Y2 = Y(J-2)
A2 = X3 - X2
M2 = (Y3 - Y2)/A2
IF (J .EQ. LO) THEN
M4 = 2*M3 - M2
GO TO 45
ENDIF
ENDIF
X5 = X(J+1)
Y5 = Y(J+1)
A4 = X5-X4
M4 = (Y5 - Y4)/A4
IF (J .EQ. 2) M2 = 2*M3-M4
45 IF (J .LE. 3) THEN
M1 = 2*M2 - M3
ELSE
A1 = X2 - X(J-3)
M1 = (Y2-Y(J-3))/A1
ENDIF
IF (J .GE. LM1) THEN
M5 = 2*M4 - M3
ELSE
A5 = X(J+2) - X5
M5 = (Y(J+2) - Y5)/A5
ENDIF
C
C --- NUMERICAL DIFFERENTIATION
C
IF (I .EQ. LP1) GO TO 52
W2 = ABS(M4-M3)
W3 = ABS(M2-M1)
SW = W2 + W3
IF (SW .EQ. 0) THEN
W2 = .5
W3 = .5
SW = 1.0
ENDIF
T3 = (W2*M2 + W3*M3)/SW
IF (I .EQ. 1) GO TO 54
52 W3 = ABS(M5-M4)
W4 = ABS(M3-M2)
SW = W3 + W4
IF (SW .EQ. 0.) THEN
W3 = .5
W4 = .5
SW = 1.0
ENDIF
T4 = (W3*M3 + W4*M4)/SW
IF (I .EQ. LP1) THEN
T3 = T4
SA = A2 + A3
T4 = .5*(M4 + M5 + A2*(A3-A2)*(M2-M3)/SA**2)
X3 = X4
Y3 = Y4
A3 = A2
M3 = M4
ENDIF
GO TO 60
54 T4 = T3
SA = A3 + A4
T3 = .5*(M1+M2+A4*(A4-A3)*(M3-M4)/SA**2)
X3 = X3 - A4
Y3 = Y3 - M2*A4
A3 = A4
M3 = M2
C
C --- DETERMINATION OF COEFFICIENTS
C
60 Q2 = (2.0*(M3-T3)+M3-T4)/A3
Q3 = (T3+T4-2.*M3)/A3**2
C
C --- COMPUTE POLYNOMIAL
C
70 DX = DX - X3
V(K) = Y3+DX*(Q1+DX*(Q2+DX*Q3))
80 CONTINUE
RETURN
END
C
C---END INTRPL
C